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Abstract —In this paper, we propose a novel design for molec¬ 
ular communication in which both the transmitter and the 
receiver have, in a 3-dimensional environment, multiple bulges (in 
RF communication this corresponds to antenna). The proposed 
system consists of a fluid medium, information molecules, a 
transmitter, and a receiver. We simulate the system with a one- 
shot signal to obtain the channel’s finite impulse response. We 
then incorporate this result within our mathematical analysis 
to determine interference. Molecular communication has a great 
need for low complexity, hence, the receiver may have incomplete 
information regarding the system and the channel state. Thus, 
for the cases of limited information set at the receiver, we pro¬ 
pose three detection algorithms, namely adaptive thresholding, 
practical zero forcing, and Genie-aided zero forcing. 

Index Terms —Molecular communication via diffusion, inter¬ 
ference, Brownian motion, 3-D simulation, symbol detection 
algorithm. 


I. Introduction 

If operating at the nano-scale is to make an impact at 
the macro scale, there must be high cooperation among 
multiple devices 0. 0- Therefore, nanonetworking emerges 
as a new paradigm and molecular communication, as an 
interdisciplinary branch, draws the attention. The literature 
has proposed various molecular communication systems, such 
as molecular communication via diffusion (MCvD), calcium 
signaling, microtubules, DNA micro-arrays, pheromone sig¬ 
naling, and bacterium-based communication |li|-(5J. In an 
MCvD, a number of micro- and nano-machines reside in a 
viscous environment and communicate through molecules that 
are emitted into the medium. Following the physical charac¬ 
teristics of the diffusion channel, these molecules propagate 
through the environment. Some of these molecules arrive at 
the receiver (i.e., hit the receiver) and form chemical bonds 
with the receptors on the surface of the receiver. The properties 
of these received molecules (e.g., concentration and/or type) 
constitute the received signal 0> §■ 

To increase MCvD performance, many enhancements are 
proposed in the literature such as incorporating inter symbol 
interference (ISI) mitigation techniques 0. §, using multi¬ 
ple molecule types |9j(, using multiple input multiple output 
(MIMO) techniques [lOj. The authors in ffO) , however, mainly 
focused on multiuser interference and paid scant attention to 
the ISI. In molecular communication, ISI is the main source 
of communication impairment and must be analyzed precisely. 
Therefore, this paper focuses on introducing and enhancing 


MIMO analysis. 

We propose a molecular MIMO system and, considering 
the channel and interference models, introduce new detection 
algorithms. The interference in the system is caused by the ISI 
and inter link interference (ILI). This is where ISI originates 
from the previous emissions of the corresponding antenna 
and ILI arises from the other antenna. To the best of our 
knowledge, this is the first work that considers MIMO for 
molecular communication while taking into account the ISI 
and ILI. We model the channel’s impulse response by mod¬ 
ifying the single input single output (SISO) channel model 
in a 3-dimensional (3-D) environment EJ- Consequently, 
we analyze the system performance via the MIMO simulator 
developed through MATLAB. 

The rest of this paper is organized as follows. In Section [II] 
we describe the system model including topology, propaga¬ 
tion, and communication models. In Section [IIIJ we detail 
the channel estimation method and the proposed detection 


algorithms. In Section IV the results and discussions are 
presented. Finally, the conclusion is given in Section [V] 


II. System Model 


Consider a molecular communication system in a 3-D 
environment with two point sources and two spherical receiver 
antennas Q The transmitter releases, without any directional 
preferences, a certain number of messenger molecules at 
once. Without colliding into one another or undergoing any 
chemical reactions, the emitted molecules travel along the fluid 
medium via diffusion. When a molecule hits the boundary of a 
spherical antenna, it is immediately absorbed by the receiving 
antenna and removed from the medium. We assume that the 
transmitter-receiver pair is synchronized and the receiver can 
count the number of received molecules during a symbol 
duration. 

During communication, the transmit antennas convey inde¬ 
pendent messages to their corresponding receiver bulges. Each 
link of transceivers uses the same type of molecule and the 
molecules from the other transmitter cause ILL 


A. Topology and Propagation Model 

As shown in Fig. [T] there are two transmitter antennas, 
Txl and Tx2, placed d distance apart from the corresponding 

'Throughout this paper, we use the terms bulge and antenna interchange¬ 
ably. 
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Fig. 1. Topological model of molecular 2x2 MIMO system. 

receiver bulges, Rxl and Rx2. We have two receiver bulges 
with the same radius r r , which are placed h distance apart. 
The centers of the receiver bulges, Txl, and Tx2 lie on the 
corners of a rectangular grid. The receiver bulges are attached 
to the receiver body and we assume that just the antennas are 
capable of receiving molecules, which is a realistic assumption 
considering many examples found in the nature. For example, 
epithelial cells, neurons, and migrating cells are examples of 
polarized cells that have heterogeneous receptor deployments, 
which is an adaptation to the environment and the signaling 
mechanism. 

When a molecule is released from a point source, its 
movement in a fluid is governed by diffusion and drift. Drift 
is applicable if there is a flow and we leave the drift case to 
future work. The dynamics of diffusion can be described by 
Brownian motion. Derived in CD is the fraction of molecules 
that are absorbed by a single spherical receiver antenna in a 
3-D environment until time t with given parameters 

d - D) = wfi erfc (jm) (1) 

where D denotes the diffusion coefficient. In our setup, how¬ 
ever, there are two receiver bulges and we cannot directly use 
0 Therefore, we simulate the Brownian motion for released 
particles within the given MIMO setup by using 

(x t , y t , z t ) = (xt-At, yt-At,z t -At) + (Ax, Ay, A z) 

Ax ~ 7V(0, 2D At) 

Ay ~ 7V(0, 2D At) (2) 

Az ~ AT(0,2D At) 

where Xt, yt, zt, and AT(p,cr 2 ) are the particles’ positions at 
each dimension at time t, and the normal distribution with 
mean // and variance cr 2 . The Brownian motion simulator for 
the MIMO setup is a modified version of the simulator that is 
developed for a SISO case in J7]|. During its trip, if a molecule 
hits one of the spherical receiver bulges, then it is absorbed 
and removed from the environment. Therefore, a molecule 
can contribute to the signal just once and we have four 
different F(t\r r ,d, D) values depending on the molecule’s 
emission source and hitting bulge. We use the simulation data 


to formulate F(t |ry, d, D) for the 2x2 molecular MIMO setup 
and utilize it for the analysis. 

B. Communication Model 

To encode information, we use the binary concentration 
shift keying (BCSK) modulation technique. We denote Q\ 
and Qq as the number of molecules released to send bit-1 
and bit-0, respectively. The transmitter has independent sets 
of bits X\ and x 2 for their own messages. During the n th 
symbol, Txl and Tx2 send x\ [n] and x 2 [n] each by releasing 
Q Xl [ n ] and Q x ,,\ n molecules at the start of the symbol time, 
and wait until the next emission time. The duration between 
consecutive symbols is called the symbol duration and is 
denoted by t s . The number of molecules absorbed at the 
receiver stochastically follows a binomial distribution with a 
hitting probability, which is related to t s , d, ry, and D |2j. We 
define F,j (L , t 2 ) as the probability of hitting to Rxi between 
t\ and t 2 for a molecule that is released from a Txj. So we 
can define the random variable Sij(t\. t, 2 ) as follows: 

Sij(ti,t 2 ) — B ( 1 , Fij(ti,t 2 )), ( 3 ) 

where B(n,p) denotes the binomial random variable with n 
trials and success probability p. Binomial random variable 
in can be considered as Bernoulli trial. S l3 is utilized 
while evaluating the number of received molecules at R x i that 
originates from Txy. 

In this paper, we consider two types of interference sources 
for the receiving bulge: interference from the previous symbols 
of the corresponding transmitter (i.e., ISI) and interference 
from the current and the previous symbols of the other link 
(i.e., ILI). The ISI at the n th symbol can be modeled as a sum 
of interference due to the previous symbols. The ILI at the n th 
symbol can be modeled as a sum of interference due to the 
other link emissions including the current symbol. Hence, the 
interference model at Rxi can be expressed as 

Ii[n] = ISIi[n] +ILIi[n] 

n—1 

ISIi\p\ ^ ^ Qxj\n—k]^ii ffi] 

k =1 
n—1 

ILI 1 [n\ = ^2 Qx 2 [n-k\Si2[k\ ( 4 ) 

k—0 

n—1 

IL^2\p\ = ^ ^ Qxi\n—k]^2l[^\ 
k—0 

Sij[k] A Sij(kt s , (k+l)t s ) 

where Ii[n], ISI t [n], and ILIi\n\ denote random variables 
of total interference, ISI, and ILI induced at R x i at the n th 
time slot respectively. Note that the summation in the ILI term 
starts from zero, because the current symbol of the other link 
also induces interference. The channel output at Rxi for the 
n th time slot can be written as 

Vsxi[n] = {(Qi-Qo)xi[n\ + Qo)Sii[0]+I z [n\+rii[n\ (5) 











where, for each time slot, yyi xi is the random variable of 
the number of received molecules and n, denotes the molec¬ 
ular noise that can occur due to outer molecular invasion 
or decomposition at Rxi. The effect of noise is assumed 
to be a normal distribution J\f(p n: cr^). In this paper, the 
transmitter emits zero molecules to send bit -0 to reduce the 
energy consumption and separate the signal amplitudes as 
much as possible. Therefore, substituting Q o in (|5]i with zero 
and writing in matrix form for each symbol slot yields the 
following 
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for 2 x 2 MIMO system. We write 0 in short as 


Fig. 3. Representation of detection algorithms and their requirements. 


y = Hx + I + n. (7) 

The details of I will be given in the following section. 

III. Fitting Channel Parameters and Proposed 
Detection Algorithms 

We need to have a formula for the first hitting probability in 
molecular MIMO setup to formulate Therefore, exten¬ 

sive simulations are carried out to understand the underlying 
formula and we use nonlinear curve fitting on the simulation 
data. After having the fitted cumulative distribution functions 
(CDFs) we carry out our analytical derivations via utilizing 
approximate F, :) (t) functions. We also introduce the proposed 
detection algorithms and optimal thresholds in this section. 

A. Fitting Channel Parameters 

We use a model function that is coherent with 0 to 
fit the simulation data (i.e., the formula is coherent with 
molecular SISO system in a 3-D environment with some 
control parameters). The model function structure for non¬ 
linear fitting is as follows: 

F — ((|r -' d ■ D) = Sy erfc ((IcFis) ,8) 

where by, b 2 , and b : > are controllable parameters. We run 
extensive simulations for different parameter sets and estimate 
mean CDF of the hitting molecules. The hitting molecules are 
separated according to where they are originated for finding 
Fn, F\ 2, F 21 , and F 22 . Due to the symmetry of the topology, 
Fyy and F 12 are very close to F 22 and F 2 \ in the simulation 
data. 

NonLinearModel class in MATLAB is used for fitting non¬ 
linear regression models. We implement the model function 
and utilize the simulation outputs to have a closed form CDF 
estimation obeying (| 8 j. Fig. [2] depicts the distance versus fitted 
model parameters for different h and ry values. The values of 
b 2 and 63 change little when the distance is increased and, 
similar to the SISO case, they are close to 0.55. 


B. Detection Algorithms 

In this section, we introduce four detection algorithms. 
Each of them requires a different set of information given 
to the receiver. Fig. [3] describes the algorithms in terms of 
the information required. There is a default set that is needed 
commonly for all the algorithms. The default set consists of 
system parameters D, t s , and topology parameters such as d, 
h, and ry. First algorithm works with only the default set and 
we name it the fixed threshold method. It uses a predetermined 
threshold and does not adapt to varying t s or Qy. The second 
algorithm, called adaptive thresholding, additionally requires 
Qy, which is not a big assumption since Qy is determined with 
the communication protocol and the modulation. The output 
of the detector, y a , is formulated in ([9]) and the algorithm 
calculates the optimal decision threshold accordingly. 

Va = = yj-l-Hz + Z + n) (9) 

VI VI 

The third algorithm, called practical zero forcing method, 
needs the average channel response matrix, which is denoted 
by H. Inspired from the zero forcing of the conventional 
communication strategy, the formulation is given in ( flO) . The 
output of the detector, y p , is formulated as 

y p =H~ 1 Hx + H~ 1 I + H~ 1 n. (10) 

Lastly, the receiver is aided by Genie so that the exact channel 
states are known for every signal reception time. We call the 
algorithm for this case Genie-aided zero forcing. The output 
of the detector, y (J , is formulated in (ITTk. It provides the best 
performance but is not feasible sincethe randomness of the 
molecular communication channel originates from the Brown¬ 
ian motion of molecules and is hard to acquire instantaneously. 

y g =x + H 1 I + H 1 n. (11) 

Our main contribution includes proving that the adaptive 
thresholding and the practical zero forcing methods perform 
exactly the same for the symmetrical MIMO topology. 

Theorem 1. When the centers of the transmitter and the 
receiver antennas form a rectangular grid, the detector outputs 
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of the adaptive thresholding and the practical zero forcing 
methods satisfy y a = A 0 y p (Aq denotes the hitting probabil¬ 
ities of a molecule to the intended bulge at the current symbol 
duration). 


approximate I\[n\ = ISIi[n\ + ILIi[n], given in (|4|, with 
a Gaussian distribution that has mean /// and variance a 2 . 
Lemma [2] provides the formulations for estimating the mean 
and the variance of the interference. 


Proof: The topological symmetry guarantees the random 
variables <Sn[fc] and <S22[fc] have equal statistical parameters 
for a positive integer k. They are both binomial random values 
with success probability of A k . The transmitter sends Q\ 
molecules for transmitting a bit-1. Therefore, the diagonal 
entries of H follow a binomial distribution and approximated 
to the normal distribution as follows: 


H(i, i ) ~ B(Qi,Ao) « A f ( Q 1 A 0 , QiAq(1 — Ao)) (12) 

and the channel mean H equals Q 1 A 0 E where E denotes a 
2x2 identity matrix. It leads H = (1/QiAq)E and y a in 
(|9j» becomes a multiplication of Ao and y p in (10 1 . ■ 

Theorem [T] ensures that both methods on average perform 
the same, since the signal detection properties (i.e., the detector 
outputs) are similar up to a constant multiple. 


C. Interference Formulations and Optimal Thresholds 

In this section, we formulate the interference and find 
the optimal decision rule for the adaptive thresholding and 
practical zero forcing methods analytically. To derive the 
optimal decision threshold, the study uses the maximum-a- 
posterior (MAP) method. 

As the receiver is unaware of the parameter Q \ for the 
case of fixed threshold method, we have to predetermine a 
range of Q\ to use and decide the decision threshold, rjf, 
to minimize the bit error rate (BER) for all Q i. Note that 
the analysis of Genie-aided zero forcing inherits difficulties of 
acquiring instantaneous H. Hence, the optimal thresholds for 
Genie-aided zero forcing are found empirically. 

We should consider the interference when determining the 
thresholds and consequently the received symbol. Topological 
symmetry ensures that all the properties of Rxl and Rx2 
coincide in terms of interference, so it is sufficient to analyze 
Rxl. The channel output for Rxl at n th time slot becomes 
?/RxiM = QiSn[A\xi[n\ + h[n\ + m[n] from (|5j. We 


Lemma 2. The k th ISI term in the sum¬ 

mation of (|?]) has mean value of TtiQ\A k and variance of 

niQiA k {l - A k ) + 7ro7riQfA%. 

Proof: With probability 7^ and 7r 0 , Q Xl [n-k]S\\[k\ fol¬ 
lows JV(QiA k ,Q 1 A k (l — A k )) and becomes just zero, re¬ 
spectively. Therefore, the mean of the received ISI becomes 
TtiQ\A k and the variance becomes 

<j 2 = E[x 2 } — E[x } 2 

= TTi(Q\A k + QiAl(l — A k ) 2 ) — w{Q 2 A k 
= 7TlQfAl(l — A k ) 2 + (iTl — TTi)QfA k 
= 7t 1 Q 1 A k (l - A k ) +7r 0 7riQ 2 A 2 . 


In a similar way, we can apply the lemma to find the mean 
and variance of the k th ILI and sum both to find the total 
mean and variance of I\[n\. The total interference mean and 
variance at the n th symbol slot becomes 


f n— 1 


n— 1 


pi = t t 1 q 1 ( A k + B k ) 

\k—l k—0 ) 

( n— 1 n— 1 

A k +^ bi 

k—l k—0 

/ n— 1 i 

+ Tt\Qi I A k ( 1 — A k ) + B k ( 1 — B k ) 


C 2 I 


n— 1 


(13) 


\k =1 


k =0 


where B k denotes the success probability of both <Si 2 [fc] and 
S 21 \k]. Note that ( fl3| l does not require the previously sent bit 
sequences. It requires only the index of the current symbol, 
as it evaluates the expected value over the cases. Hence, © 
can be used for each symbol consecutively. 

After formulating the interference and the detector output, 
we can now derive the thresholds for the practical zero forcing. 





























We denote the probability density function of y p when the 
transmitted bit is 0 and 1 as y p | 0 and y p | l5 respectively. 
We approximate the detector outputs at Rxi by Gaussian 
distribution as 


Vp |o(*) o,ctq) 

V P \i(i) ~ -AT(/^i, crj) 


= J\f 
= J\f 


( 

WiAi ’ 
^1+po, 


+ <fj 

<9Mo 
(i ~ At) 
QiA 0 



(14) 


where y p \ 0 (i) is z th element of the 2x1 vector. Note that y p ^ 
is the sum of y p | 0 and H H that determines the mean and 
the variance of y p \\(i) in ( fl4| . The formulations are obtained 
by utilizing ( fl0l > and ( | 12 [ i. 

Now we define a decision rule as argmax(y p | i ) and need 
to find intersection points of two distributions (i.e., to find the 
decision threshold, rj p , for practical zero forcing ). It leads to 
the equality 


1 ( (Vp~do) 2 \ _ 1 

W2^ XP V 2<Tq ) ai V^ 

and its solution in terms of y p becomes 


(Vp-lLl ) 2 

2a\ 



Vp = Mo + 


-1± Vl + Qg-lXl + qgfilnfl) 

/3-1 


for p = (cti/cto ) 2 > 1. We denote the bigger one as 77 + and 
the smaller one as r]~, then the decision rule for the decoded 
bits x becomes: 


* = S(y p ) 


0 Vp > y'v > v P 

1 otherwise 


where 8(.) is the decision function at the receiver. 

Note that /3 > 1 because af is the sum of erg and the 
variance of H H. The case where 3 becomes 1 means that 
y p | 0 and y p ^ have the same variances and it is trivial that the 
threshold becomes Bo±Bk_ 

We can obtain similarly the decision rule for adaptive 
thresholding. We define y a | 0 and y a ^ as in and their 
means and variances, by Theorem [l] become A 0 p 0 , (A 0 <to) 2 > 
A 0 p 1 , and (Aqct^ 2 , respectively. Therefore, we have an 
equation similar to © for finding the decision threshold of 
adaptive thresholding method 

A 0 cri f {vw-A 0 po) 2 \ ( {Vw — Aopi) 2 \ 

^ 0^0 6XP { 2 (A 0 <7o) 2 )- C ^{ 2(^ 0( ti) 2 ) 

that can be solved similarly. 


IV. Results 

The system parameters used in this paper are given in 
Table [I] We first give the definition of signal-to-interference- 
ratio (SIR) metric in molecular MIMO system and analyze 
the effect of distance, ry, and h. Next, we use the BER as the 
performance metric and analyze the effect of Qi and t s . 


TABLE I 

Range of Parameters Used in the Analysis 


Parameter 

Variable 

Values 

Diffusion cefficient 

D 

50 iirn z / s 

Distance 

d 

(2, 4 }fim 

Radius of the receiver 

7> 

(2, 4 }fim 

Bulge separation 

h 

{1, 2 }/im 

# molecules for sending 
bit-1 

Qi 

{100 rsj 600} molecules 

Probability of sending 
bit-1 

1X1 

0.5 

Symbol duration 

t s 

{0.05 ~ 1} sec 

Molecular noise variance 


10 

Bit sequence length 


5 x 10 4 

Replication 


20 



Fig. 4. SIR plots of different topologies (D = 50 pnr 1 /s) 


A. SIR Analysis 

SIR is defined as the ratio of the expected number of 
molecules coming from the intended transmitter in the in¬ 
tended time slot and mean ILI plus ISI for just a one-shot 
signal. 

_ _ -^11 (0, t a ) _ 

-Fii(f s ,oo) + Fi 2 (0,oo)' 

Note that this definition is specific to the molecular communi¬ 
cation case and explains the clearness of the mean signal term 
in the received signal. 

Fig. [4] depicts t s versus SIR for different topology parame¬ 
ters. With respect to SIR, the best enhancement is provided by 
reducing the distance between the transmitter and the receiver. 
Increasing the bulge size also gives merit, however, using 
more separation distance between Rx bulges results in only 
insignificant improvement. Thus, we can conclude that SIR 
is influenced more by the ISI term than the ILI term for our 
system setup. 

We set the topological parameters d = 2pm, r> = 4 pm, 
h = 2 pm, and D = 50pm 2 /s for the rest of the performance 
evaluation. The selected system parameters and the fitted 
values for model parameters are given in Table [II] Utilizing 
fitted values enables us to estimate Fij(t) analytically. 

















































TABLE II 

Fitted model parameters for the selected topology (d = 2 fim, 
r r = 4 fim, h = 2fim, D = 50 


Function 

b i 


f>3 

Tn(t) 

0.9155 

0.5236 

0.5476 

F 12 {t) 

0.1534 

0.2780 

0.5363 



Fig. 5. BER performance of detection algorithms (f„ = 80 ms). 


B. BER Analysis 

In this section, we analyze the BER with respect to varying 
Q i and t s and compare the performance gain of each of the 
proposed detection algorithms. Each Tx sends 5 x 10 4 bits 
with equal probability of sending bit-1 and bit-0. Most prior 
work has shown that, with an appropriate symbol duration, the 
current symbol is mostly affected by one previous symbol with 
the rest being negligible @. Therefore, we consider 

four slots of interference in the simulations. We examined 
all the thresholds between —1 and 2 with KT 3 interval and 
checked the optimal threshold for each Q\. Empirically found 
fixed threshold, rjf, is selected as 0.2. 

Fig. 0 shows the BER performance of detection algorithms 
while Q i varies from 100 to 600. The first observation is 
that the adaptive thresholding and the practical zero forcing 
provide coinciding results, as in Theorem [T] Increasing Q\ 
(i.e., the signal power) decreases the BER for all the detection 
algorithms. However, the improvement of fixed thresholding is 
significantly lower than the other methods. When the instan¬ 
taneous H is known, Genie-aided zero forcing is applicable 
and gives the best performance. Obtaining that information, 
however is not easy. On the other hand, obtaining the optimal 
threshold by knowing Q i and 7iq is feasible and leads to a 
performance that is close to Genie-aided zero forcing. 

Fig. [6] illustrates the BER performance against t s from 
50 ms to 130 ms. It shows that increasing t s gives faster 
improvement in terms of BER relative to increasing Q\. This 
means that to achieve a lower BER it is more effective to 
decrease the information rate than to increase the transmit 
power. 



Fig. 6. BER performance of detection algorithms (Q i = 500). 


V. Conclusions 

Data rates in molecular communications are acutely affected 
by interference. Therefore, any enhancement for molecular 
communications should consider interference effects precisely. 
In this paper, we proposed a MIMO system for the MCvD 
that takes into account inter symbol and inter link interference. 
Moreover, we proposed three symbol detection algorithms that 
depend on the information set that the receiver has. First, we 
modeled the channel’s finite impulse response via fitting 3- 
D MIMO simulator results considering the fraction of the 
received molecules. We utilized the estimated function (that 
gives the fraction of received molecules) to determine inter¬ 
ference and the optimal thresholds for the proposed methods. 
In the performance analysis, we first analyzed the effect of 
varying topological conditions on the SIR. The result shows 
that the transmitter-receiver distance and the size of receiver 
bulges (antennas) are more effective at reducing the interfer¬ 
ence rather than the separation of bulges. We investigated the 
performance of the proposed detection algorithms in terms of 
BER while varying Q i and t s . We quantified the enhancement 
of these parameters in the molecular MIMO system. As a 
future direction, we will focus on a testbed implementation 
of the proposed algorithms while incorporating the drift. 
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